function [H, Kernel] = CreateKernel(beta, mu)

if(0)

  Nx = 1024; 
  lc = 16;
  dx = 1 / Nx;
  x = (0:Nx-1)'/Nx;

  epsilon = 1/Nx;

  lambda = 1d-2;
  lambdarnd = 1d-2;
  lambdamod = 3d-1;




  % Purely random
  if(0)
    V = lambdarnd * rand( Nx, 1 );
  end

  % Filtered random noise
  if(0)
    V = lambdarnd * rand( Nx, 1 );
    fV = fft(V);
    kcutoff = 30;
    ind = find(abs(k)>kcutoff);
    fV(ind)=0;
    V = real(ifft(fV));
    V = V / max(abs(V)) * lambda;
    %  V = V * Nx / kcutoff;
  end

  % Periodic 
  if(0)
    V = zeros(Nx,1);
    % V = lambda * cos( 2*pi / lc * (x/epsilon) );
    % mu = 1.904722371438797e-02;
    mu = 1.921471959676978e-02;
  end

  % Periodic + random
  if(0)
    V = lambda * cos( 2*pi / lc * (x/epsilon) ) + lambdarnd * rand( Nx, 1 );
  end

  % Periodic + simple long wave mode
  if (0)
    V = lambda * cos( 2*pi / lc * (x/epsilon) )  + ...
      lambdarnd * cos(2*pi * x );
  end

  % modulated periodic
  if(0)
    V = lambda * cos( 2*pi / lc * (x/epsilon) ) .* (1 + ...
      lambdamod * cos(2*pi * x ) + lambdamod * sin(2*pi *2*x)); % +  0.01 ...
    %      * randn(size(x)));
  end

  e = ones(Nx,1);
  H = spdiags(0.5*[-e 2*e -e], -1:1, Nx, Nx) / dx^2  * epsilon^2; % Discretization of the laplacian
  H(1,Nx) = -0.5/dx^2*epsilon^2;
  H(Nx,1) = -0.5/dx^2*epsilon^2;

  H = H + spdiags(V,0,Nx, Nx);
  % D = eigs(H, 100, 'sm' );

  Kernel = beta * (H - mu * speye( Nx ) ) + 1i * pi * speye(Nx);
end

Nx = 1024; 
dx = 1 / Nx;
x = (0:Nx-1)'/Nx;

lambda = 1d-2;
lambdarnd = 1d-2;
lambdamod = 3d-1;

  if(1)
    V = lambdarnd * rand( Nx, 1 );
    fV = fft(V);
    kcutoff = 80;
    ind = find(abs(fV)>kcutoff);
    fV(ind)=0;
    V = real(ifft(fV))*10;
  end


e = ones(Nx,1);
H = spdiags([-e 2*e -e], -1:1, Nx, Nx) / dx^2; % Discretization of the laplacian
H(1,Nx) = -1/dx^2;
H(Nx,1) = -1/dx^2;
H = H + 1 * spdiags(ones(Nx,1),0,Nx,Nx) + spdiags(V,0,Nx,Nx);
Kernel = H;
